December 16, 2022
Features
Overview
Toy example
Actual case study: The Mitchell River, Australia
Shiny Application
Reproducibility, Academic use, …
Source: https://medium.com
Powered by Rcpp and RcppArmadillo
Find optimal solutions
Solve a specific new problem using mathematical programming
Ideally, fast 🙂
Three questions: Planning objective, Optimization objective, and Intensity of threats.
Bolch, E.A. et al. (2020). Remote Detection of Invasive Alien Species
Source: https://ib.bioninja.com.au
inputData())Data that defines the spatial prioritization problem (planning units data, feature data, threats data, and their spatial distributions).
inputData( pu, features, dist_features, threats, dist_threats, sensitivity = NULL, boundary = NULL )
What if I have my data in another format?
A: There are functions that allow you to import and convert data to data.frame format
Examples:
##from .txt to data.frame
df <- read.table('myfile.txt',sep='\t')
##from .csv to data.frame
df <- read.csv("myfile.csv")
Source: https://libguides.umn.edu
pu)data.frame
Specifies the planning units (PU) of the corresponding instance and their corresponding monitoring cost and status. Each row corresponds to a different planning unit. This file is inherited from the pu.dat in Marxan.
id
integer unique identifier for each planning unit.
monitoring_cost
numeric cost of including each planning unit in the reserve system.
status (optional)
integer value that indicate if each planning unit should be available to be selected (0), locked-in (2) as part of the solution, or locked-out (3) and excluded from the solution.
## id monitoring_cost status ## 1 2 0 ## 2 2 0 ## 3 2 0 ## 4 2 0 ## 5 2 0 ## 6 2 0
pu)Game, E. T. y H. S. Grantham. (2008). Manual del Usuario de Marxan: Para la versión Marxan 1.8.10
pu)Source: https://marxansolutions.org
pu)0)2)3)Source: https://marxansolutions.org
features)data.frame
Specifies the conservation features to consider in the optimization problem. Each row corresponds to a different feature. This file is inherited from the spec.dat in Marxan.
id
integer unique identifier for each conservation feature.
target_recovery
numeric amount of recovery target to achieve for each conservation feature.
target_conservation (optional)
numeric amount of conservation target to achieve for each conservation feature.
name (optional)
character name for each conservation feature.
##id target_recovery name ## 1 11 feature1 ## 2 16 feature2 ## 3 8 feature3 ## 4 9 feature4
features)features)Source: https://ecology.fnal.gov
features)features)features)dist_features)data.frame
Specifies the spatial distribution of conservation features across planning units. Each row corresponds to a combination of planning unit and feature.
pu
integer id of a planning unit where the conservation feature listed on the same row occurs.
feature
integer id of each conservation feature.
amount
numeric amount of the feature in the planning unit. Set to 1 to work with presence/absence.
#pu feature amount #1 3 1 #2 3 1 #3 3 1 #4 3 1 #5 3 1 #6 3 1
dist_features)Source: https://calumma.typepad.com
dist_features)Source: https://ib.bioninja.com.au
threats)data.frame
Specifies the threats to consider in the optimization exercise. Each row corresponds to a different threats.
id
integer unique identifier for each threat.
blm_actions (optional)
numeric penalty of connectivity between threats. Default is 0.
name (optional)
character `name for each conservation feature.
#>id name blm_actions #> 1 threat1 0 #> 2 threat2 0
threats)Source: https://docs.marxanweb.org
dist_threats)data.frame
specifies the spatial distribution of threats across planning units. Each row corresponds to a combination of planning unit and threat.
pu
integer id of a planning unit where the threat listed on the same row occurs.
threat
integer id of each conservation feature.
amount
numeric amount of the threat in the planning unit. Set to 1 to work with presence/absence. Continuous amount values require that feature sensitivities to threats be established.
#>pu threat amount action_cost status #> 8 2 1 2 0 #> 9 2 1 2 0 #>10 2 1 2 0 #>11 1 1 3 0 #>11 2 1 4 0 #>12 1 1 3 0
dist_threats)data.frame
specifies the spatial distribution of threats across planning units. Each row corresponds to a combination of planning unit and threat.
action_cost
numeric cost of an action to abate the threat in each planning unit.
status (optional)
integer value that indicate if each planning unit should be available to be selected (0), locked-in (2) as part of the solution, or locked-out (3) and excluded from the solution. :::
dist_threats)Source: https://ib.bioninja.com.au
sensitivity)data.frame
specifies the sensitivity of each feature to each threat. Each row corresponds to a combination of feature and threat.
feature
integer id of each conservation feature.
threat
integer id of each threat.
a (optional)
numeric the minimum intensity of the threat at which the features probability of persistence starts to decline.
b (optional)
numeric the value of intensity of the threat over which the feature has a probability of persistence of 0.
#>feature threat #> 1 1 #> 2 1 #> 3 1 #> 4 1 #> 1 2 #> 2 2
sensitivity)data.frame
specifies the sensitivity of each feature to each threat. Each row corresponds to a combination of feature and threat.
c (optional)
numeric minimum probability of persistence of a features when a threat reaches its maximum intensity value.
d (optional)
numeric maximum probability of persistence of a features in absence threats. :::
sensitivity)sensitivity)boundary)data.frame
Specifies the spatial relationship between pair of planning units. Each row corresponds to a combination of planning unit.
id1
integer id of each planning unit.
id2
integer id of each planning unit.
boundary
numeric penalty applied in the objective function when only one of the planning units is present in the solution.
#>id1 id2 boundary #> 1 1 0 #> 2 1 1 #> 3 1 2 #> 4 1 3 #> 5 1 4 #> 6 1 5
minimizeCosts) vs Maximize Benefits (maximizeBenefits)\(x_{ik}\) is the decisions variable that specifies whether an action to abate threat \(k\) in planning unit \(i\) has been selected \((1)\) or not \((0)\).
\(c_{ik}\) is the cost of the action to abate the threat \(k\) in the planning unit \(i\).
\(b_{is}(x)\) is the benefit of the feature \(s\) in the planning unit \(i\) (ranging between \(0\) and \(1\)).
\(t_s\) is the recovery target for feature \(s\).
\(f(x)\) is the function of connectivity penalty.
minimizeCosts) vs Maximize Benefits (maximizeBenefits)\(x_{ik}\) is the decisions variable that specifies whether an action to abate threat \(k\) in planning unit \(i\) has been selected \((1)\) or not \((0)\).
\(c_{ik}\) is the cost of the action to abate the threat \(k\) in the planning unit \(i\).
\(b_{is}(x)\) is the benefit of the feature \(s\) in the planning unit \(i\) (ranging between \(0\) and \(1\)).
\(B\) is the budget available.
\(f(x)\) is the function of connectivity penalty.
problem())Create an optimization model for the multi-action conservation planning problem, following the mathematical formulations used in Salgado-Rojas et al. (2020).
problem( x, model_type = "minimizeCosts", budget = 0, blm = 0 )
solve())Solves the optimization model associated with the multi-action conservation planning problem. This function is used to solve the mathematical model created by the problem() function.
solve( a, solver = "", gap_limit = 0, time_limit = .Machine$integer.max, solution_limit = FALSE, cores = 2, verbose = TRUE, name_output_file = "output", output_file = TRUE )
solve())gurobi, cplex, symphony)Rsymphony
solve())Source: https://www.math.uwaterloo.ca
solve())getActions()Returns the spatial deployment of the actions for each planning unit of the corresponding solution.
getActions(x, format = "wide")
Output example,
## solution_name pu 1 2 conservation connectivity ## 1 sol 1 0 0 0 0 ## 2 sol 2 0 0 0 0 ## 3 sol 3 0 0 0 0 ## 4 sol 4 0 0 0 0 ## 5 sol 5 0 0 0 0 ## 6 sol 6 0 0 0 0
getSolutionBenefit()Returns the total benefit induced by the corresponding solution. The total benefit is computed as the sum of the benefits obtained, for all features, across all the units in the planning area.
getSolutionBenefit(x, type = "total")
Output example,
## solution_name feature benefit.conservation benefit.recovery benefit.total ## 1 sol 1 0 11 11 ## 2 sol 2 0 16 16 ## 3 sol 3 0 10 10 ## 4 sol 4 0 9 9
getCost()Provides the sum of costs to actions and monitoring applied in a solution.
getCost(x)
Output example,
## solution_name monitoring threat_1 threat_2 ## 1 sol 61 20 65
Features
Toy example
Overview
Actual case study: The Mitchell River, Australia
Shiny Application
Package prioriactions can be found at CRAN, so it can be installed using:
install.packages("prioriactions")
Latest stable versions can be downloaded and installed from GitHub as follows (package remotes should be installed first):
if (!require(remotes)) install.packages("remotes")
remotes::install_github("prioriactions/prioriactions")
So, we load the prioriactions package.
# load package library(prioriactions)
prioriactions website)This example contains 100 planning units, 4 features and 2 threats.
The distribution of features and threats can be plotted on a grid of 10 x 10 units.
prioriactions contains this simulated example inside the setup files. You can extract it by the data() function.
The monitoring cost values ranging from 1 to 10 and all status of 0 (not locked).
# load planning unit data from prioriactions data(sim_pu_data) #To load simulated data #head(sim_pu_data)
A RasterLayer object can be used to present this spatial information. The pixel values correspond to the monitoring costs of each planning unit.
library(raster) #To plot rasters r <- raster(ncol=10, nrow=10, xmn=0, xmx=10, ymn=0, ymx=10) values(r) <- sim_pu_data$monitoring_cost plot(r)
Contains information about the features such as its id and targets (mandatory when minimizing costs).
# load features data from prioriactions data(sim_features_data) #head(sim_features_data)
Contains information on the spatial distribution of these features across planning units.
# load features data from prioriactions data(sim_dist_features_data) #head(sim_features_data)
To plot the spatial distribution of the first feature,
# load amount of features data
features <- reshape2::dcast(sim_dist_features_data,
pu~feature,
value.var = "amount",
fill = 0)
values(r) <- features$`1`
plot(r)
Provides information about the threats such as their id and name.
# load threats data from prioriactions data(sim_threats_data)
Provides information on the spatial distribution of these threats
# load threats data from prioriactions data(sim_dist_threats_data)
Indicates which features is sensitive to what threat.
# load threats data from prioriactions data(sim_sensitivity_data)
Provides information on the spatial relationship between PU and they are presented in long format.
# load boundary data from prioriactions data(sim_boundary_data)
After having loaded our data, we will now create the data object through the inputData() function.
# create conservation problem
b <- inputData(pu = sim_pu_data,
features = sim_features_data,
dist_features = sim_dist_features_data,
threats = sim_threats_data,
dist_threats = sim_dist_threats_data,
sensitivity = sim_sensitivity_data,
boundary = sim_boundary_data)
# print problem
print(b)
## Data ## planning units: data.frame (100 units) ## monitoring costs: min: 1, max: 10 ## features: feature1, feature2, feature3, feature4 (4 features) ## threats: threat1, threat2 (2 threats) ## action costs: min: 1, max: 10
It is advisable to determine the maximum benefit achievable for each conservation feature through the getPotentialBenefit() function
# get benefit information getPotentialBenefit(b)
# create optimization problem c <- problem(b, model_type = "minimizeCosts")
## Warning: The blm argument was set to 0, so the boundary data has no effect
## Warning: Some blm_actions argument were set to 0, so the boundary data has no ## effect for these cases
# print model print(c)
## Optimization Problem ## model sense: minimization ## dimensions: 284, 396, 12.656 kB (nrow, ncol, size) ## variables: 396
To solve our model, we need an optimizer. In this case, we will use the Rsymphony library.
# solve optimization problem d <- solve(c, solver = "symphony", verbose = TRUE, output_file = FALSE, cores = 2) # print solution print(d)
## Solution overview ## name: sol ## objective value: 738 ## gap: 0% ## status: Optimal solution (according to gap tolerance: 0) ## runtime: 0.017 sec
getActions()# get actions from solution actions <- getActions(d, format = "wide") # print first six rows of data # head(actions)
getActions()values(r) <- actions$`1` plot(r)
# create optimization problem c2 <- problem(b, model_type = "minimizeCosts", blm = 10)
## Warning: Some blm_actions argument were set to 0, so the boundary data has no ## effect for these cases
# print model print(c2)
## Optimization Problem ## model sense: minimization ## dimensions: 29984, 10296, 883.856 kB (nrow, ncol, size) ## variables: 10296
# solve optimization problem d2 <- solve(c2, solver = "symphony", verbose = TRUE, output_file = FALSE, cores = 2) # print solution print(d2)
## Solution overview ## name: sol ## objective value: 863 ## gap: 0% ## status: Optimal solution (according to gap tolerance: 0) ## runtime: 22.847 sec
# get actions from solution actions2 <- getActions(d2, format = "wide") values(r) <- actions2$`1` plot(r)
Features
Overview
Toy example
Actual case study: The Mitchell River, Australia
Shiny Application
We divided the whole catchment (71,630 km2) into 2316 sites (i.e., sub-catchments)
We sourced the distribution of 45 fish species in the Mitchell river catchment as our conservation features.
We considered four major threats to freshwater fish species in the catchment: water buffalo (Bubalis bubalis), cane toad (Bufo marinus), river flow alteration (caused by infrastructure for water extractions and levee banks) and grazing land use.
Features
Overview
Toy example
Actual case study: The Mitchell River, Australia
Shiny Application
Source: https://shiny.rstudio.com